research papers 



Received XXXXXXX 0000 
Accepted XXXXXXX 0000 



o 
o 

(N 
>> 



0^ 
(N 



> 

in 
o 
o 



I 

o 
o 



New phasing method based on the principle of 
minimum charge 

Pavel Kalugin^* 

'^Laboratoire de Physique des Solides, Universite Paris-Sud, 91405 Orsay France. Correspondence e-mail: 
kalugin@lps.u-psud.fr 

A new method of the phase determination in X-ray crystallography is proposed. 
The method is based on the so-called "minimum charge" principle, recently sug- 
gested by Elser The electron density function p is sought in the form p(x) = 
|'0(x)p, where ip is an «-component real function. The norm / \ijj{x)\^dx is min- 
imized under the constraint imposed by the measured data on the amplitudes of 
Fourier harmonics of p. Compared to the straightforward implementation of the 
"minimum charge" scheme, the method attenuates the Gibbs phenomenon and is 
also capable of extrapolation of the diffraction data beyond the set of measured 
amplitudes. The method is applicable to quasicrystals under the condition that the 
number of components n of the function ip is bigger than the dimensionality of 
the "atomic surface". It is successfully tested on synthetic data for Fibonacci chain 
and the octagonal tiling. In the latter case the reconstructed density map shows the 
shape of the atomic surface, despite relatively low data resolution. 
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1. Introduction 



The vast majority of the existing direct methods of X-ray struc- 
ture determination approach the phase problem as a problem of 
constrained minimization. The quantity to minimize plays the 
role of the likelihood functional, optimization of which is sub- 
ject to constraints imposed by the known structure factor am- 
plitudes. The choice of this quantity is a matter of tradeoff be- 
tween three requirements. First, this functional should approx- 
imatively represent the common notion about the likelihood of 
a given density function. Second, it should be effectively com- 
putable, and finally it should allow for efficient minimization. 
Usually none of these requirements is fully satisfied. Consider 
for instance the traditional direct methods based on the condi- 
tional probability distribution of the structure invariants (Harker 
& Kasper, 1948; Karle & Hauptman, 1950). The derivation of 
the expressions for the conditional probability takes as a start- 
ing point the assumption that unconditional probability of the 
density distribution (its Bayesian prior) is a translationally in- 
variant measure on the ensemble of (5-like atoms in a unit cell 
(Cochran, 1955). Clearly this approach misses the physical con- 
straint on a minimal distance between the atoms. Although this 
prior measure allows one to express the conditional probability 
of the phase invariants in a closed form (Cochran, 1952; Haupt- 
man & Karle, 1953; Cochran, 1955), the resulting formulae are 
effectively computable only for invariants of small order. As a 
result, instead of the true conditional probability function the 
practical algorithms use some sort of ad hoc approximation, 
usually based on combination of the triplets and quartets (see 
e.g. "Shake 'n Bake" algorithm (Weeks et al., 1993)). However, 
even these simplified functionals are not easy to optimize. The 
existing phase refinement methods are based on iterative proce- 
dures and can easily get stuck in a local minimum of a func- 
tional. 

Recently, Elser (1999) suggested the so-called principle of 



minimum charge. According to this principle, the correct set 
of phases should minimize the average charge density (the un- 
known Fourier component with k = 0, which has to be added to 
the density to make it non-negative). The rationale behind this 
principle is that the functions satisfying it tend to have shallow 
highly degenerated minima and spiky maxima, which is what 
one would expect of an atomistic density function. Compared to 
the solid statistical foundations of the conventional direct meth- 
ods, the principle of minimum charge may look arbitrary. How- 
ever, this comparison is not fair, because on the way from the 
first principles to the practical implementation of the conven- 
tional methods many ad hoc assumptions are made. In contrast, 
the principle of minimum charge is almost readily applicable in 
the phase refinement algorithms. Additional advantage of this 
principle is that it could be used without modification for the de- 
termination of structure of quasicrystals and non-commensurate 
crystals. This is especially important, because the conventional 
methods fail when applied to these structures. Indeed, a naive 
attempt to approximate quasicrystals by crystals with very big 
number of atoms in the unit cell gives rise to the divergence of 
the normalized structure factors En — 0{N^^^) and of the triplet 
amplitudes Ahk = 0{N). This immediately leads to meaning- 
less results (e.g. that phases of all triplets Ohkl = 0). 

A straightforward implementation of the principle of mini- 
mum charge (e.g. that proposed by Elser (1999)) implies solv- 
ing a minimax problem. Indeed, the set of phases <I>k satisfies 
the principle if the deepest minimum of the "density with zero 
average" 

p'(r) =2V|FK|cos(K-r-<I>K) (1) 



K 



over r has the maximal possible value over all possible sets Ox- 
Robust algorithms for finding a global saddle point are much 
more difficult to design than those finding a global minimum 
or maximum. The reason is that the methods which are usually 
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applied to prevent the algorithm from getting stuck with a lo- 
cal minimum (e.g. simulated annealing or multiple runs) do not 
work in the case of local saddle points. In this article we pro- 
pose a version of the principle of minimum charge which can 
be formulated as a problem of global minimization, and not as a 
minimax problem, allowing one to circumvent the shortcomings 
of a straightforward approach. 

The naive implementation of the principle of minimum 
charge suffers from another drawback, which is due to the so- 
called Gibbs phenomenon. In order to understand this problem, 
consider the case when the true values of all phases <E>k are 
known. Suppose also, that the correction for the atomic form 
factors is included in the values of the structure factors Fk in 
(1), or in other words that the atoms are point-like. Due to the 
Gibbs phenomenon, any finite sum of the form (1) presents neg- 
ative "bumps" around the positions of atoms. The depth of these 
artifacts may be as big as 20% of the peak bight. This picture is 
clearly different from multiply degenerate shallow minima that 
one would expect of a function satisfying the principle of min- 
imum charge. In other words, the correct set of phases gives a 
density function (1) which is suboptimal from the point of view 
of this principle. Any further optimization of it may only intro- 
duce errors in the values of phases. The proposed method sig- 
nificantly attenuates the importance of the Gibbs phenomenon, 
although does not remove it completely. 

2. Method 

2.1. Representation of the density function 

The cornerstone of the new method is the representation of 
the atomic density function. The traditional way of reconstruct- 
ing this function consists of using finite Fourier sums (1) with 
possible application of weights aimed to attenuate Gibbs phe- 
nomenon. A different approach is used here. The density func- 
tion is approximated by a square of a finite Fourier sum. More 
precisely, the density function p(x) is modeled as 

p(x) = |v(x)|^ (2) 

where V'(x) is a multi-component real function. Each compo- 
nent ipa of ip has a finite Fourier spectrum: 

^^„(x) = ^^„,Ke'''■^ (3) 
KeA 

where A is a finite subset of reciprocal lattice vectors. The com- 
plex coefficients V'a.K satisfy the condition 

V'a.-K = i>a,K (4) 

Before going any further, let us discuss the rationale behind 
this representation. Recall, that the density function satisfying 
the principle of minimum charge should possess the following 
properties: 

1. The amplitudes of its Fourier components for the set 
of measured Bragg reflections should correspond to the 
measured data ^k. 



2. The global minimum of the function should be highly de- 
generate (see Fig. 1). 
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Figure 1 

Phase refinement based on the principle of minimum charge leads to 
the reconstructed density maps with highly degenerate shallow min- 
ima. 

The traditional way to represent the density function by a finite 
Fourier sum automatically guarantees that the condition 1 is sat- 
isfied. However, the condition 2 cannot be expressed in a closed 
form as an explicit constraint on the phases <t>K. On the other 
hand, degenerate global minimum occurs naturally in the func- 
tions of the type (2) when ip{x) has multiple zeros. In the same 
time, the condition 1 can be explicitly imposed on the function 

(2) . This constraint can de formulated as a system of equations 
of 4-th order on the variables 'ti>a,-K (see Appendix A). 

Mention should be made of the criterion for choice of the sup- 
port A for the Fourier spectrum of the function ip in the equation 

(3) . Clearly, if the set A is too small, the condition 1 may be im- 
possible to satisfy. This happens e.g. if there are wave vectors 
L belonging to set E of measured data, which cannot be repre- 
sented as L = H — K, where H, K e A. On the other hand, 
choosing the set A equal to the set S guarantees that any con- 
straint on the amplitudes of Fourier components of density pK 




with K e S can be satisfied in the representation (2). Indeed, 
assuming the foUowing constraint on the amplitudes of pK- 

\Pk\ = OK, 

setting 

f ifK = 0, Q=0 

V'a.K = < foK ifK^O, a = 

[ if a 7^0 

gives the correct Fourier components of p in the limit e — > -|-0: 

\pk\ =aK + 0{e). 
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Note, that in this case the support of the spectrum of p in this 
case is roughly twice as large in the reciprocal space as the set 
A. In other words, the representation (2) allows for extrapola- 
tion of the structure factors. One can further extrapolate the ex- 
perimental data by choosing an even larger set A. The limit of 
possible improvement of the resolution in this way is set by the 
occurrence of spurious peaks in the density map in the case of 
over-extrapolation due to the increased sensitivity of the result 
to the errors in the measured data. 



2.2. Optimization criterion and constraints 

The representation (2) of the density is a non-negative func- 
tion. Hence the average value of p can be used as a figure of 
merit for the optimization, in other words, the principle of min- 
imum charge regains its original meaning. It should be empha- 
sized, that instead of the problem of finding a global saddle 
point (global maximum of global minima) we deal here with an 
easier problem of constrained minimization. The role of param- 
eters is played by the Fourier components VS^.k of the function 
tf}. The average density is expressed through ^/S^^k as 



KeA a 



(6) 



The parameters VSq-.k can be considered as a real vector in the 
Mn-dimensional space, where M is the number of points in the 
set A and n is the number of components of '0- The principle of 
minimal charge thus boils down to the problem of constrained 
minimization of the norm of this vector 

Consider now the constraints imposed by the measured data 
on the values V'a,K- Let us assume that the amplitudes of the 
structure factors |Fk| are measured for the wave vectors K be- 
longing to a subset S of the reciprocal lattice. Then the seem- 
ingly obvious constraint on the density function p(x) from (2) 
consists of the requirement that 



for any K e S. 



This condition expressed in terms of the Fourier components of 
tj) from (3) takes the form: 



a H,K-HeA 



|Fk| foranyKeS. (7) 



This constraint, however, does not take into account the effect 
of the Gibbs phenomenon. The importance of this effect is clear 
from the Figure 2. 




Figure 2 

Artifacts occurring due to Gibbs phenomenon. The original structure 
is a one-dimensional crystal with one atom per unit cell. The number 
of independent structure factors is 10. 

This plot shows the function p(x) obtained by the minimiza- 
tion of the average density (6) under the constraints (7). The 
structure factors used here correspond to a one-dimension crys- 
tal with one atom per unit cell, that is |Fk| = 1 for all K e S. 
Clearly the phasing on the Figure 2 is incorrect. The reason is 
that the correct set of phases would correspond to a function 
with deep negative bumps around the "atom" because of trun- 
cation of data in the reciprocal space. The optimization routine 
attempts to flatten these bumps out and align all minima on the 
same level, giving rise to a wrong solution. 

The standard approach to reduction of the Gibbs phenomenon 
consists of using the windowing in the Fourier domain. As ap- 
plied to the formula (7) the windowing consists of multiplying 
the ampUtudes of the structure factors by weights: 



E 

a H,K-HeA 



(8) 



The choice of the coefficients wk is a matter of tradeoff between 
softening the truncation of data on the reciprocal lattice and pre- 
serving the spatial resolution of the density map. In the field of 
digital signal processing there exist many windowing formulae 
designed to reduce the magnitude of Gibbs phenomenon, e.g. 
Welch or Hann windows (Rabiner & Gold, 1975; Bloomfield, 
1976). Nonetheless, instead of using these formulae, we shall 
derive the expressions for wk which are optimized from the 
point of view of the discussed method. Ideally, the weights wk 
should guarantee that the global minimum of the average den- 
sity (6) with the constraints (8) corresponds to the density map p 
with the correct phases. The weights should be a function of the 
sets S and A only, and should not depend on the values of the 
structure factors Fk. Unfortunately, such ideal weights do not 
exist. This is clear at least from the fact that the perfect phasing 
is impossible when the set of measured data is too small as com- 
pared to the number of atoms in the unit cell. Nevertheless, as 
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long as the resolution of the measured data is sufficiently high, 
there exists a nearly optimal set of weights. More specifically, 
the weights wk chosen in such a ways as to guarantee the cor- 
rect phasing for the structure with one atom per unit cell should 
produce reasonable results for other structures as well, provided 
that the peaks in the density map do not overlap. 

Let us start the construction of the optimal weights with 
the case of a one-dimensional crystal with one atom in a unit 
cell. First of all, note that the constraints in the formula (7) 
with |Fk| = 1 are equivalent to those of the equation (8) with 
l^kl = 1 /wk- In other words, one can interpret the density map 
on the Fig. 2 as a result of correct phasing of an a priori un- 
known structure factor amplitudes \Fk\ — 1/wk- Clearly, the 
number of peaks in this density function depends on the number 
of constraints in (7). One would naturally expect that the mini- 
mization of the global charge with only one constraint leads to 
some very simple structure. Indeed, as shown below, if the set 
E in (7) contains only one wavevector Kq, which is equal to the 
elementary period of the reciprocal lattice, the minimization of 
the average density gives rise to a structure with one atom in a 
unit cell. Let 'ipa,K = Va.K be the solution of (7) with the single 
constraint \Fk^ | = 1 • One can use these values to construct the 
following set of weights: 



EkgaEo l'?a,K| 



(9) 



expression 



+ y E (12) 



k=-N 



k=-N+l 



Differentiation with respect to 'ip gives the equation for the point 
of extremum: 



(13) 



(here we assume 'tp-N-i = '4>n+\ = 0). In other words, the solu- 
tion of the principle of minimum charge with the constraint (10) 
should be an eigenstate of a lattice Laplace operator with zero 
boundary conditions beyond the set A. Straightforward calcula- 
tions show that the following function 



ijjk=C cos 



2{N+ 1) 



(14) 



This set of weights guarantees that the density map obtained 
with a single constraint will satisfy the principle of minimum 

1/2 

charge for any number of constraints. Indeed, V'a.K = fja,K 
obeys the equations (8) with the weights (9) and |Fk| = 1 for 
all K G E. In other words, the weights wk are optimal for the 
set A. 

Let us actually compute the optimal weights for the case 
when the set A includes the nodes —A', — + 1 , ... A' of the 
one-dimensional reciprocal lattice. We suppose that ip{x) has 
only one component and use simplified notations by writing -tpi^ 
or F), instead of 'tjja,K or Fk (here k is the number of the node). 
In the case when only the amplitude of the structure factor F\ is 
known, there is only one constraint on the values {i^,}: 



-N+l 



= wi\Fi\ 



(10) 



By introducing Lagrange multiplier the constrained minimiza- 
tion of the average density can be replaced by finding an uncon- 
strained extremum of the following quantity: 



Yl '^ki>*k + A 



k=-N 



k=-N+l 



(11) 



Note, that the second sum in the formula (11) can be always 
made real and positive by an appropriate shift 6x of the coordi- 
nate system and multiplication of ipk by e'*^''^. Then if expression 
(11) has an extremum at a given the same is true for the 



gives the smallest value of the average density. The density is a 
sum of peaks centered at the lattice nodes 



pix)=J2f[i^N + 2){x-n) 



(15) 



nez 



where the shape of an individual peak is given by the following 
formula: 



f[u] = C 



/ cos(7rM) 
V 1 - 4m2 



(16) 



As one might expect, in the limit A/^ — > oo the density map tends 

to a sum of (^-functions corresponding to a crystal with one atom 
per unit cell. The formula (9) gives the values for the optimal 
weights: 



sin {(2N + 3 - \k\)0) + (2N + 3 - |/v|) sm((9) cos{ke) 

sin {{2N + 3)0) + {2N + 3) sin(6l) cos(A;6l) ' 

(17) 

where ^ = 7r/2(A'^-|-l). The results of testing of these weights 
with the structure factors corresponding to one atom in the unit 
cell are shown on Figure 3. 
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Figure 3 

Reconstructed density map for the same test case as that of the Fig. 2, 
but with constraints (8). The values of weights are given by the formula 

(17) . 

Consider now the ways to generalize the above construction 
of the optimal weights to more than one dimension. In this case 
the minimal number of constraints needed to define a structure 
is equal to the dimensionality d of the crystal. Solving the prob- 
lem of minimum charge with the known structure factors for 
K = Ki . . . Krf yields to the formula similar to (13): 

V^a,K + y (V'a,K-K, + V^a,K+K,) = (18) 

i 

^a,K = forK^A. 

The weights can be obtained from the solution •il>a,K of this 
equation using the formula (9). 

Contrary to the one-dimensional case, equation (18) does not 
yield to a unique set of optimal weights. First of all, the choice 
of the wave vectors K, is ambiguous. By analogy with the one- 
dimensional crystal, it seems natural to require that the vectors 
K, form a basis of the reciprocal lattice. However, there are 
many ways to choose a basis of a lattice in more than one di- 
mension. The other problem is related with the choice of La- 
grange multipliers A,. In the formula (13) the solution ip^ is al- 
ways one of the eigenstates of Laplace operator and the role of 
the parameter A is restricted to mere selection of the eigenstate. 
In contrast, V'a.K in (18) may vary continuously with A,. In or- 
der to resolve the ambiguity one needs to recall that the density 
corresponding to the solution of the finite difference equation 

(18) describes the structure with one atom at the origin of the 
unit cell. In other words, the function V'(x) should have a sharp 
peak at origin, that is its Fourier components 'tpa,K should vary 
slowly across the domain A. The requirement that the equation 
(18) should admit of such slowly varying solution determines 
the choice of parameters K; and A,-. 



2.3. Quasicrystals and incommensurate structures 

Consider now the apphcation of the discussed method to the 
deterministic quasicrystals and incommensurate structures in 
general. These structures can be conveniently described using 
the so-called superspace or "cut-and-project" method (de Wolff 
etal., 1981; Duneau & Katz, 1985). According to this approach, 
the density function of quasicrystals can be obtained as a cut 
through a periodic function in a space of higher dimensionality. 
The Fourier spectrum of the structure is obtained as a projection 
of the spectrum of the periodic function on the cut direction, and 
thus consists in a discrete sum of (5-functions. If the direction of 
the cut is incommensurate with the periodicity, each ^-peak in 
the Fourier spectrum of the quasicrystal corresponds to a node 
of the reciprocal lattice of the periodic function. By this means 
the phase problem for quasicrystals can be reformulated in a 
conventional way as a problem of phasing for a periodic func- 
tion in a space of higher dimensionality. 

For the case of point-like atoms the density function of a qua- 
sicrystal in the real space is a discrete sum of ^-functions. The 
corresponding periodic density function consists of S-like dis- 
tributions on sub-manifolds, commonly referred to as "atomic 
surfaces" (Janssen, 1986; Bak, 1986). This brings up again the 
question of optimality of the weight factors wk from (8). In- 
deed, the weights obtained following the method described in 
the previous section are appropriate for distributions of narrow 
non-overlapping peaks. Atomic surfaces clearly do not fall into 
this category. As a result one could expect that using the weight 
factors optimized for point-like atoms might give rise to incor- 
rect results in the case of atomic surfaces. Nevertheless, lacking 
any better alternative, we applied these weights for quasicrys- 
tals as well. Despite concerns, the results of the numerical tests 
described below show no significant distortion of the density 
map. 

Another problem emerges when the atomic surfaces are dis- 
continuous, which is the case for all known deterministic struc- 
ture models of real materials. The point is that due to the finite 
resolution in the reciprocal space, the boundary of the atomic 
surface is inevitably smoothed out. As a result, when the phys- 
ical space cuts through an atomic surface near its boundary, the 
amplitude of the corresponding peaks in the density map is re- 
duced. Similarly, when the cut misses the atomic surface by a 
short distance, a "phantom" peak appears in the density map. In 
reality both of the above phenomena occur simultaneously, be- 
cause of the so-called "closedness" property of the atomic sur- 
faces (Katz, 1989). This results in occurrence of double peaks of 
reduced intensity, which are often separated by a distance much 
smaller than the typical interatomic spacing. A careful analy- 
sis shows that such double peaks in places group together to 
form more complex clusters (Kalugin & Katz, 1993). It should 
be emphasized that all these artifacts are not specific to the 
phasing method and will exist in any reconstructed quasicrys- 
talline density map. They should not be confused with true par- 
tially occupied sites which occur in real quasicrystals because 
of the phason disorder (Lyonnard et ah, 1996; de Araujo et ah, 
1996; de Boissieu et al, 1995). 
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2.4. Choice of the number of components of 

To this point the number of components of the function il>{x) 
has been of no importance for us. However, because of the iter- 
ative nature of the optimization algorithm, this parameter plays 
an important role in the case of quasicrystalline structures. The 
problems which arise in the case when the number of compo- 
nents of tp is too small are clear from Fig. 4. This figure de- 
picts a one-dimensional atomic surface in the two-dimensional 
space crossed by a line of zeros of one-component function tp. 
The resulting hole can only be removed by pushing it towards 
the boundary of the atomic surface. It may occur, however, that 
pushing the hole inwards decreases the average density. In this 
case, the algorithm will converge to a local minimum. 




Figure 4 

An unremovable hole in the atomic surface on the reconstructed den- 
sity map. Such holes occur when a Une of zeros of one-component 
function ip crosses the atomic surface. 

The above problem could be avoided if the function tp had 
two components. In this case, the zeros of tjj are points. Would 
such point fall onto the atomic surface, it could be easily pushed 
away by a small perturbation of tp. The same applies to the case 
of more than two dimensions. In this case the number of com- 
ponents of tp should be bigger, than the dimensionality of the 
atomic surface. 

3. Results 

The results presented in this section were obtained numerically 
using the algorithm of constrained minimization described in 
the Appendix A. The algorithm was tested on the sets of syn- 
thetic data including one-dimensional crystal and quasicrystal, 
as well as a two-dimensional quasicrystalhne structure. No ro- 
tational symmetry was assumed in any of the cases. When the 
structure actually possessed additional symmetry, it was recov- 
ered as a result of optimization. 

Because of its iterative nature, the algorithm of Appendix A 
does not guarantee convergence to the global minimum of the 
average density. It is worth noting, however, that for all tested 
structures the correct phasing corresponds to the deepest of the 
found minima. 

3.1. One-dimensional crystals 

The test structure includes five point-like atoms. Their 
'charges' and fractional coordinates x are given in the Table 1 . 



As the structure does not possess central symmetry, the recon- 
structed density can correspond to any of the enantiomorphs 
with equal probabihty. 
Table 1 

One-dimensional test structure 
'charge' x 
2.0 0.0 
1.5 0.6 
1.2 0.25 
1.0 0.43 

LO 08 

The density map shown on the Fig. 5 has been obtained with 
the first 12 structure factor amplitudes. The function ip had 2 
components. The support of the Fourier spectrum of each com- 
ponent A from the formula (3) included the wavevectors from 
the interval [—12, 12]. The positions of five peaks on the Fig. 5 
correspond to the coordinates of atoms in the Table 1 up to 
a global translation. The amplitudes of peaks are also qualita- 
tively recovered. One can remark, however, that the two small- 
est peaks on the Fig. 5 are noticeably different, although they 
correspond to identical atoms. The value of the average den- 
sity is equal to 6.636, which is about 1% smaller than the total 
charge 6.7 of the Table 1. The iterative optimization converged 
to the global minimum solution in 14 of 20 trials when start- 
ing with random normally distributed 'ipa,K was 70%, giving the 
70% success rate. 




Figure 5 

The restored density map corresponding to the structure from Table 1 . 
The input data includes first 12 structure factor amplitudes. The func- 
tion ip has two components and is defined with the same resolution. 

The test structure from the Table 1 was also used to bench- 
mark the extrapolating capabilities of the algorithm. The input 
data was restricted to the first 9 structure factor amplitudes. 
Note, that the structure recovery with any smaller data set would 
be impossible because a one-dimensional structure of 5 atoms 
is described by 9 real parameters not including a global trans- 
lation. By contrast, the set A from (3) included all wavevectors 
from the interval [—50,50]. By this means, the algorithm ex- 
trapolated the "experimental" scattering data to a region of the 



research papers 



reciprocal space roughly 1 1 times larger that that where they 
were initially defined. The reconstructed density map is shown 
on the Fig. 6. The positions of atoms are perfectly accurate up 
to a global translation and inversion. Note also, that despite a 
smaller resolution of the input data (which is only about 0.65 of 
the smallest interatomic distance), the amplitudes of the peaks 
are restored with higher accuracy than on the Fig.5. The value of 
the average density is equal to 6.690, which is also much closer 
to the total charge 6.7 of the Table 1. This can be explained 
by overlapping of wider peaks on the Fig.5, which makes the 
correct solution suboptimal. As a result, the algorithm produces 
slightly distorted solutions, which have smaller values of the 
average density. Narrower peaks on the Fig.6 make this phe- 
nomenon much less visible. 
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Figure 6 

Density map for the test structure of the Table 1 extrapolated from the 
low resolution data. 



One could anticipate a lower success rate in the latter exam- 
ple because of the bigger number of parameters in ^c«,k to fit. 
Contrary to the expectations, convergence to the global mini- 
mum was obtained in 19 trials of 20. 

3.2. One-dimensional Fibonacci chain 

Fibonacci chain is probably the best known example of a one- 
dimensional quasicrystal. The chain consists of identical atoms 
arranged on a straight line. The distance between neighbouring 
atoms may take two values, usually referred to as "short" and 
"long" segments of the chain. The ratio of lengths of short and 
long segments equals r = (\/5 — l)/2. The short and long seg- 
ments alternate following a quasiperiodic Fibonacci sequence, 
which is usually defined recursively (Senechal, 1995). The same 
structure can be obtained following the conventional "cut-and- 
project" method (see Fig. 7). 




Figure 7 

The "cut-and-projcct" method of construction of the Fibonacci chain. 
The one-dimensional physical space cuts the two-dimensional lattice of 
periodically arranged segments ("atomic surfaces"). The atoms, shown 
as filled circles, are located at intersections of "atomic surfaces" with 
the physical space. 

The algorithm has been tested on a set of 40 independent re- 
flections with the Miller indices hP' + Jc'' < 25. Although the 
structure on the Fig. 7 possesses central symmetry, this symme- 
try was not imposed as an additional constraint on tp. The recon- 
structed density map is shown on the Fig. 8. The optimization 
has yielded the correct solution with 100% success rate for 50 
trials with random initial conditions. 




Figure 8 

Reconstruction of the atomic surfaces of the one-dimensional Fi- 
bonacci chain. The input data includes 40 independent reflections with 
+ < 25. The function tj; has two components. The shown area 
of the density map includes 4 unit cells. The axes represent fractional 
coordinates. 
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3.3. Two-dimensional octagonal Ammann-Beenker 
tiling 

Ammann-Beenker octagonal tiling (Ammann et ai, 1992; 
Beenker, 1982) refers to quasiperiodic covering of a plane by 
squares and 45° rhombi. Structure models associated with this 
tiling can be constructed by decorating each tile with atoms. 
The simplest decoration consists of placing atoms at vertices of 
squares and rhombi. The resulting structure can be obtained by 
"cut-and-project" technique from a periodic density function in 
four-dimensional space. The corresponding atomic surfaces are 
two-dimensional perfect octagons, one per unit cell. 

The input data included the reflections with h^+k^+l^+m^ < 
5 (where h, k, I and m stand for four-dimensional Miller in- 
dices). There are 68 pairs of opposite nodes of the reciprocal lat- 
tice satisfying the above inequality. The 8mm symmetry of the 
diffraction pattern makes only 14 of them independent. How- 
ever, as mentioned above, the point symmetry was not taken 
into account. As a result all 68 pairs of reflections were con- 
sidered as independent. The function ^ had three components. 
Note, that as no central symmetry is imposed, the coefficients 
V'q.k from the formula (3) are complex numbers. 

The optimal set of coefficients ^pa.K returned by the algo- 
rithm should be converted to the density map in the physical 
space. As the direction of the physical space is incommensurate 
with the four-dimensional periodic lattice, this conversion im- 
plies Fourier transform of unevenly spaced data. This precludes 
using standard FFT algorithms, which significantly slows down 
displaying the results. A workaround for this problem consists 
of tilting the physical space slightly to make it commensurate 
with the lattice. This is equivalent to replacing the quasicrystal 
by a close approximant (Duneau & Audier, 1994). 
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Figure 9 

Density map corresponding to the vertices of an approximant to 
Ammann-Beenker tiling, reconstructed from the reflections with h' + 

^ ^ ^- < 5. 



Figure 9 depicts the reconstructed density map for an ap- 
proximant to Ammann-Beenker tiling. The approximant is ob- 
tained by replacing the vectors Ui = (\/2/2, 1/2,0,-1/2) and 
U2 = (0, 1 /2, V2/2, 1 /2) spanning the physical space by the 
rational vectors: 



u'l - (5/7,1/2,0,-1/2) 
= (0,1/2,5/7, 1/2). 



(19) 



Note that the peaks on the density map substantially overlap 
because of low resolution of the input data. This is also con- 
firmed by the fact that the algorithm yields the average den- 
sity as low as 0.718 of its correct value, suggesting significant 
over-optimization. Nevertheless, most of the atomic positions 
are correctly resolved, as can be seen from comparison with the 
vertices of the ideal approximant tiling shown on Fig. 10. The 
double peaks which are expected due to the smoothened edges 
of the atomic surface (see the explanation above) are not fully 
resolved. They manifest themselves as elongated features on the 
density map. The optimization algorithm has converged to the 
correct solution with 100% success rate of 10 trials with random 
initial conditions. 



k' + l' 



Figure 10 

The vertices of the approximant to Amman-Beenker tiling with the 
physical space spanned by the vectors (19). 



One could also consider cutting the four-dimensional den- 
sity function along different directions. For instance, a cut in 
the direction of the atomic surface would reveal its shape. In 
practical conditions this shape will be distorted because the cut 
may not pass exactly through the center of the four-dimensional 
peak corresponding to the atomic surface. Such cut is shown on 
Fig. 1 1 . Despite low resolution of the input data one can clearly 
see a faceted octagonal shape. 
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Figure 11 

The cut through the reconstructed density map of Ammann-Beenker 
tiling in the direction parallel to the atomic surface. The scale is arbi- 
trary. 

4. Summary and discussion 

We have described a new method of structure determination 
based on the minimum charge principle. The method possesses 
extrapolating capabilities, and can restore atomic positions from 
low resolution data (about 65% of the interatomic distance). It is 
also applicable to quasicrystals and incommensurate structures. 
As the current implementation of the algorithm was intended 
as a demonstration of principle, no special efforts were taken 
to make it more efficient. Making the new method practical for 
crystals with non-trivial point symmetry and for real quasicrys- 
tal requires further research. 

The failure to take advantage of the symmetry of the structure 
is a major drawback of the current implementation. The prob- 
lems related to the symmetry are the most visible in the case of 
central symmetric structures. It is well known, that in the case 
of central symmetry, the structure factors Fk are real numbers 
(the same appUes to some structure factors when the point sym- 
metry group contains other elements of order 2). As a result, 
each of the constraints (8) on the absolute values of the Fourier 
components of p{x) defines two disconnected manifolds in the 
nM-dimensional space of coefficients -tpa.K- The total number 
of disconnected manifolds defined by all constraints is equal to 
2™, where m is the number of real structure factors. In its cur- 
rent implementation, the optimization algorithm sticks to one of 
these pieces in the early stages and and then continues looking 
for the point of minimal average density on this piece only. In 
this way, the global minimum will most likely be missed. 

The other problem is related to the fact that the symmetry of 
the function tp{^) does not necessarily coincide with that of the 
density function pix). In the general case, rotations and trans- 
lations in the real space are accompanied by orthogonal trans- 



formations in the «-dimensional space of components of tp. In 
other words, the symmetry group of ip{x) is an extension of 
the symmetry group of p{x) by a subgroup of 0{n). Generally 
speaking, the search for the solution corresponding to the global 
minimum of the average density should be performed over all 
such extensions. 

Appendix A. 
Numerical algorithm 

This section describes in details the algorithm used to min- 
imize IV'P with the constraints (8). The algorithm converges 
quadratically in the vicinity of a local minimum with the com- 
putational complexity of one iteration 0{{Mn)^) and the stor- 
age requirement of the order 0{{Mn)^). Here M stands for the 
number of points in the set A in (3) and n is the number of 
components of tp{x). These properties are suitable for refining 
an already found approximate minimum. In other words, this 
algorithm should be used as a second stage in a two stage op- 
timization scheme. Nonetheless, the results of testing described 
above show that this algorithm can work also in a single stage, 
starting with random values of parameters. 

The minimal charge problem can be stated in the following 
abbreviated form: 

minimize 

subj ect to /j,- ( V') = Ci . (20) 

Here ip is a vector the real Mn-dimensional space. The mini- 
mization of its norm is subject to nonlinear equality con- 
straints (20) representing the conditions (8). Note, that by squar- 
ing the both sides of (8) /z, can be chosen in the form of uni- 
form polynomials of 4-th degree in i/j. By this means all deriva- 
tives of hi are readily available in an analytic form. This makes 
appropriate application of the sequential quadratic program- 
ming (SQP) techniques of optimization (Biggs, 1975; Gill et al, 
1984). We have used the following scheme for an elementary 
SQP iteration: 

1. Compute the SVD decomposition of the Jacobian J = 
V^/i of constraints: J = USV, where UU^ = 1, VV'^ = 
1 and 5 is an X nM diagonal matrix. Denote the first 

rows of V by Vrange and the resting rows by Vnuu. 
These matrices give the projectors correspondingly onto 
the range and null spaces of the constraints. 

2. Compute Newton step in the range space: 

<5V'range = V4ge5-lf/^-(c-/j), 

where stands for a square N x N diagonal matrix of 
inverse singular values. 

3. Compute the row vector of approximate Lagrange multi- 
phers: 
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4. Compute the Hessian of the Lagrangian function L = 
ip^ip + Xh: 

dipdip 

5. Compute the projection of H onto the null space of the 
constraints //nuu = VnanHYj^^. 

6. Find the search step Stpnuw = Vj^^^x in the null space. Ide- 
ally, this step should minimize the quadratic approxima- 
tion to the variation of the Lagrangian: 

Note that the quadratic form Hnuu is not necessarily posi- 
tive definite. 

7. Update tp: 

:= V + ^V'nuU + ^V'range 

8. Repeat the steps 1-7 until convergence. 

The storage requirements of the algorithm are dominated by the 
necessity to store a dense nM x nM Hessian matrix H and its 
projection for the steps 4-6. The speed bottleneck are the steps 
5 and 6, both requiring (9((nM)^) multiplications. The SVD de- 
composition at the step 1 and the computation of Hessian at the 
step 4 take 0{N^{nM)) and 0{N{nM)'^) multiphcations corre- 
spondingly. Taking into account that M scales as A^, the con- 
tribution of these steps to the computation time may be non- 
negligible for small n. 

The only non-trivial point in the algorithm consists of han- 
dhng possibly non-positive definite Hessian matrix at the step 6. 
The Hessian often has non-positive eigenvalues when the trial 
point is far from the local minimum. Note also, that for the con- 
sidered problem the local minimum is always degenerate be- 
cause of translational invariance in the real space and rotational 
invariance in the space of the components of the field iJj{x). As 
a result, the Hessian at the local minimum always has at least 
d + n—l zero eigenvalues. The standard approach to this prob- 
lem consists of modifying the Hessian matrix to make it posi- 
tive definite. We use the fact that diagonahzation of a symmetric 
matrix is numerically stable to represent //nuu in the form 

i 

where Vi form an orthogonal basis in the null space and /x, are 
the eigenvalues of i/nuii- The modified Hessian matrix is then 
given by the formula 

J^nuu = ^rBm{\ijLi\,e]vivJ , 

i 

As the matrix H^^w is dimensionless, the regularization parame- 
ter e can be set to some constant numeric value, which should be 
much bigger than the machine epsilon but still small enough to 
keep the convergence quadratic in the vicinity of the minimum. 



The matrix H^aW is then used to compute the full Newton update 
step in the null space: 

The author is grateful to E. Tatarinova for fruitful discussions. 
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